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Abstract. We study a non-Hermitian "PT— symmetric generalization of an A'^- 
particle, two-mode Bose-Hubbard system, modeling for example a Bose-Einstein 
condensate in a double well potential coupled to a continuum via a sink in one of 
the wells and a source in the other. The effect of the interplay between the particle 
interaction and the non-Hermiticity on characteristic features of the spectrum is 
analyzed drawing special attention to the occurrence and unfolding of exceptional 
points (EPs). We find that for vanishing particle interaction there are only two 
EPs of order + 1 which under perturbation unfold either into [{N + l)/2] 
eigenvalue pairs (and in case of Af -|- 1 odd, into an additional zero-eigenvalue) 
or into eigenvalue triplets (third-order eigenvalue rings) and {N + 1) mod 3 single 
eigenvalues, depending on the direction of the perturbation in parameter space. 
This behavior is described analytically using perturbational techniques. More 
general EP unfoldings into eigenvalue rings up to (TV -|- l)th order are indicated. 
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1. Introduction 

Physical models usually describe only a rather small separated system uncoupled from 
the rest of the world. In quantum physics the behavior of such a system is governed 
by a Hermitian Hamiltonian operator. If one wants to take the coupling to some 
external world into account, one ends up with the description of an open quantum 
system. A somehow crude but instructive way to describe such open quantum systems 
is the use of effective non-Hermitian Hamiltonians. These descriptions in general 
yield complex eigenvalues whose imaginary parts describe the rates with which an 
eigenstate decays to the external world. Most often non-Hermitian Hamiltonians are 
introduced heuristically, although this approximative description can be achieved in 
a mathematically satisfactory way for example by applying the Feshbach projection 
operator technique [1] . 

Thinking of the description of open quantum systems it might be surprising 
that there is a whole class of non-Hermitian Hamiltonians which in some parameter 
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regions give rise to purely real eigenvalues and to unitary dynamics. These so-called 

T'T— symmetric Hamiltonians [2 7] possess space-time-reflection symmetry, e.g., they 
commute with the VT operator, where the operators V and T are defined by their 
effects on the position and momentum operator x and p as 

V : X —X, p I— > —p 

T : X X, p —p, i I— > — i. (1) 

In some parameter region, the region of unbroken PT— symmetry, all eigenvalues of 
"PT— symmetric Hamiltonians are purely real and the behavior of the system is similar 
to that of Hermitian quantum systems. To get a feeling for the underlying reasons 
of this "pseudo-closed" behavior of an open system, one can think of it in terms of a 
balanced probability- flow [8] . Replacing the condition of Hermiticity by the condition 
of PT— symmetry therefore yields as well a fully consistent quantum theory, which 
attracted a lot attention in the last years and actually stimulated the research in other 
fields of physics like complexified classical systems, supersymmetry and quantum field 
theory [9]. 

At first glance one might get the impression that the non- Hermiticity is nothing 
but a small perturbation which does not change the behavior of a system too much 
compared to the Hermitian case, despite an additional decay behavior, or it may 
even be equivalent to a Hermitian theory in the presence of PT— symmetry. But 
actually non-Hermitian physics can differ radically from Hermitian physics, especially 
in the presence of eigenvalue degeneracies. While a Hermitian operator is always 
diagonalizable (eigenvalues may coalesce, nevertheless they always correspond to 
distinct eigenvectors) , for a non-Hermitian Hamiltonian the occurrence of nontrivial 
Jordan blocks in its spectral decomposition is possible - there may be points 
in parameter space at which both eigenvalues and eigenvectors coalesce, so-called 
exceptional points (EPs) [10, 11]. The occurrence of EPs in a system has drastic 
effects on the systems behavior, especially concerning adiabatic features and geometric 
phases. For the occurrence of EPs in various physical models see, for example and not 
aiming at any completeness [1, 12-28]. In the theory of PT— symmetric quantum 
systems EPs naturally occur as phase-transition points between sectors of exact 
PT— symmetry and sectors of spontaneously broken PT— symmetry. 

The field of PT— symmetric Hamiltonians is still young and the underlying 
mathematical structures are not completely understood yet. Therefore in the last few 
years the interest in comparatively simple systems with a finite dimensional Hilbert 
space, especially PT— symmetric matrix Hamiltonians, was rapidly growing [23,24,29]. 
In this context most investigations focused on the mathematical behavior of simple 
matrix models, without demanding them to represent a physical system. Nevertheless 
under special conditions there are a lot of physical systems which indeed justify 
the description via a finite matrix model. In the present paper we introduce a 
PT— symmetric generalization of a prominent Hermitian matrix model, a two-mode 
Bose-Hubbard Hamiltonian, which in the case of particles acts on an iV -|- 1 
dimensional Hilbert space. 

The Bose-Hubbard Hamiltonian is a simple description of interacting bosons on 
a lattice, which only takes one state per lattice site into account. Originally the 
Hubbard model is a basic model of solid state physics, where it is mostly addressed 
in its fermionic version to describe the behavior of electrons in solids. In the last few 
years it is enjoying a renaissance in the context of Bose- Einstein condensates (BEG) in 
optical potentials. Due to the extremely low temperatures and the precise periodicity 



VT— symmetric Bose-Hubbard model 



3 



of the optical potential these systems provide the possibility of a clean experimental 
realization of many kinds of theoretical models for interacting many-particle systems. 
One prominent example is the superfluid to Mott insulator phase transition which was 
realized in a BEG in a three-dimensional optical potential [30] . 

Investigating large M-mode, A^-particle Bose-Hubbard systems quickly goes 
beyond the scope of numerical manageability. Therefore, due to its simple structure, 
the two-mode case, which one may think of describing a BEG in a double well trap, 
became a standard model [31-36]. In the present paper we introduce an effective 
non-Hermiticity to this two-mode Bose-Hubbard Hamiltonian in a "PT— symmetric 
way, which one can imagine as an additional source and sink of equal strength. A 
closely related - and slightly more physical - model would include only a sink and 
would yield complex eigenvalues with a negative imaginary part, describing a decay 
of particles. First theoretical results for this non-Hcrmitian two-mode Bose-Hubbard 
system were presented in [37]. As a possible realization one can think of a BEG 
in a double well trap, where the condensate could escape from one of the traps via 
tunneling. Another possibility would be the outcoupling of atoms from one of the 
traps via radiofrequency [38]. 

In the present paper we analyze the spectrum of the T'T— symmetric two-mode 
Bose-Hubbard system where we draw special attention to the occurrence and the 
unfolding of EPs. Numerical results are presented to illustrate the characteristic 
behavior of this unfolding. Furthermore we use perturbative methods which allow 
for analytic descriptions. Basic tools are the Le Verrier-Faddeev method [39] for 
the derivation of the coefficients of characteristic polynomials of matrices and the 
Newton-polygon technique for the extraction of the dominant powers of polynomial 
perturbations [11,40,41]. 

In detail the paper is organized as follows: In section [2] we introduce the two- 
mode Bose-Hubbard Hamiltonian and its PT— symmetric generalization and review 
the basic vocabulary according to Hermitian and non-Hermitian degeneracies. We 
discuss the analytically solvable limit of vanishing interaction in section [3] before we 
present numerical results on the spectrum for non- vanishing interaction in section [?1 
Finally, in section [S] we investigate some of the features of the spectrum previously 
found in the numerical studies analytically using perturbative methods. 



2. Bose-Hubbard model and basic non-Hermitian vocabulary 

The physical setup under consideration is a BEG in a double well potential which at 
low temperatures can be analyzed in a two-mode approximation. The corresponding 
Hamiltonian is that of a second quantized many particle system of Bose-Hubbard type 

H = e {a\ai — a^a^^ + v [a\a2 + a^a^ + - [a\ai — a\a-^ , (2) 

where a„, al are bosonic particle annihilation and creation operators for the ?— th 
mode, 2e is the on-site energy difference, v controls the single particle tunneling and c 
the interaction strength between the particles. In order to simplify the discussion we 
assume here that both v, and c are positiv^. The Hamiltonian commutes with the 

^) Note that the energy spectrum stays the same if the sign of v is altered, while it is turned upside 
down Eji — > —En if the sign of c is altered, which does not change the subsequent discussions in 
principle. Experimentally both negative and positive values of the interaction are possible and can 
actually be modulated via a Fcshbach resonance [42] . 
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N = a\ai 



(3) 



so that the total number N of particles is conserved. 

It is convenient to introduce angular momentum operators according to the 
Schwinger representation 



Ly — 
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a\ai — a\a2 



(4) 



which obey the su{2) commutation relation 

[Lx, Ly] — iLz, (5) 

and its cyclic permutations. In terms of these operators the Hamiltonian ^ assumes 
the form 

H = 2eL,, + 2vLx + 2cLl. (6) 

Thus for e,v,c € E it is an element of the universal enveloping algebrgQ W(su(2)) of 
the su{2) Lie algebra in its angular momentum I = N/2 representation. In addition 
we will often use the Lie algebra elements L± ~ Lx± iLy with commutation relations 
[L,,L±]=±L±, [L_,L+] = -2i,. 

In the standard basis of the angular momentum algebra \l,m), which can be 
defined by the relations 



L±\l,m) = T m){l ± m + l)\l,m ± 1) , 



L z\l , m) — m\l ^ m) (7) 

with / — N/2, the Hamiltonian H takes the form of a tridiagonal {N + 1) x {N + 
1)— matrix 



H 



( di + ci 

Vl-l 






2e|m|, 



Vl-l 

-1 + Cl- 






2cm 



Vl-l 

A<m<l 



+ C1-1 



Vl-l 

-di +ci J 



vVW 



(8) 



+ - m) = . 

In the following, for the non-Hermitian generalization of the Bose-Hubbard 
Hamiltonian and the structure analysis of the characteristic polynomials (see equation 
(|44p and below) another representation of the angular momentum basis in terms 
of mononomials in a complex variable will turn out to be most convenient. The 
representation ([7]) can be described in a standard way as monomials in a variable 
^ e C as (see, e.g., [46,47]) 

^l+m 

fmiO = ,„ .... = e ^i, -/ < m < / (9) 



^(/-TO)!(/ + m)! 



^) For universal enveloping algebras see, e.g., [43-45]. 
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with the normalization condition 

{l,j\l,n,)^^^^±^J J^^I^I^Ld'C^S,^, d'^:=dilmOdiReO. (10) 

Here, Hi denotes the space of polynomials in ^ of degree less or equal 21 + 1 [46,47], 
in which the SU{2) group representation acts. In the representation ([9|) the angular 
momentum operators act as first-order differential operators 

L,=(d^-l, L+ = -ed^ + 2lt i-=9« (11) 

and yield again the relations 

We note that ((71 [HI [TOl [TT|) is the standard complex irreducible representation 
(irrep) of the real Lie algebra su(2) (and the corresponding compact, simply connected 
group SU{2)) for fixed angular momentum I [46, 47]. The advantage of this 
complex {21 + 1)— dimensional irrep is its straightforward linear extendability to 
the complexification of su{2), i.e., to sl{2,C) (see proposition 4.6 in [48]). Such 
a complexification is necessarily encountered when one passes from real coefficients 
£, V, g in the Hamiltonian ^ to complex ones — as in our case when we pass 
from the Hermitian H to the non-Hermitian T'T— symmetric H by assuming s to be 
purely imaginary. Under the specific embedding su{2) ^ sl{2, C) the irrep dimension 
2/ + 1 remains fixed. For completeness, we further note that the complexification 
of the non-compact, not simply connected real SU{1, 1) yields another embedding 
sm(1, 1) ^ sZ(2,C) with corresponding extensions of the infinite dimensional su(l, 1) 
irreps [47,48]. Below we will consider boosts within the {21 + 1)— dimensional irrep of 
the su{2) ^ s^(2,C) embedding, not involving infinite dimensional sm(1, 1)— related 
irreps (and corresponding matrices of countably infinite order), i.e. we keep within the 
su{2) ^ sl{2, C) induced irrep although boosts are naturally connected with SU{1, 1) 
transformations . 

The previous considerations justify the expansion of a non-Hermitian 
generalization of the Hamiltonian ^ in the same basis ([7]) which yields a non- 
Hermitian matrix representation as a generalization of the matrix ([8]). 

In the present paper we investigate a situation, where we assume the on-site 
energy difference e of the two-mode Bose-Hubbard model ([2]) to be complex, while the 
parameters c and v are kept real. In particular, we focus on the case of a BEC in a 
symmetric double well, where the real parts of the energies of both modes are equal 
and therefore e is purely imaginary e = —i"f, 7 € M 

H = -2ijL^ + 2vL^ + 2cLI . (12) 

Physically such an imaginary on-site energy difference can be achieved by coupling 
the modes to a continuum so that they will be unstable — decaying and amplifying 
in a balanced way. Although e G iM. spoils the Hermiticity of the Hamiltonian iJ in a 
usual Euclidian Hilbert space, it nevertheless leaves H Hermitian in a Hilbert space 
with an indefinite inner product structure, i.e., in a so-called Krein space [49-56]. This 
is easily seen from the explicit matrix structure, which is essentially equivalent to ([8]) 
with the substitution e = —ij. The matrix H is not only symmetric, H = , rather 
it holds also 



H = VH^V 



(13) 
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where V is the standard involutory permutation (sip) matrix 



V 



/ 




V 1 



1 \ 






/ 



V=I 



(14) 



which is similar to an indefinite diagonal matrix 



p 2n X 2ri 



3 V 



In 






-In 



(2n+l)x(2n+l) ^ j, 











(15) 



and 



Obviously, V can be interpreted as parity operator which interchanges the a\a\- 
0302— related modes in ([2]) as 

V : \l,m) ^ r\l,m) = \l,-m). (16) 

Denoting, as usual, the involution operator of the complex conjugation ~ the time 
reversal operator of quantum mechanics - by T, where — I, and taking into 
account that 



ijt = TH^T = THT 
wc find that the Hamiltonian H for e G is T'T— symmetric 
[VT, H] = 0. 



(17) 



(18) 

According to ()13|1 it is self-adjoint in the Krein space IC-p with the indefinite inner 
product [., .Jp — (.jT'l.). Therefore the spectrum of H will contain not only real 
branches, but also pairwise complex conjugate branches and exceptional points 
(branch points) at the transitions between real and complex sectors of the spectrum. 
This is the typical behavior of a PT— symmetric operator. For completeness, we note 
that purely real branches correspond to parameter regions of exact PT— symmetry 
(H and its eigenfunctions are T^T— symmetric), whereas pairwise complex conjugate 
eigenvalues correspond to regions of spontaneously broken PT— symmetry (in contrast 
to H its eigenfunctions are not PT— symmetric). 

In the more general case of complex on-site energy difference e G C where neither 
real nor imaginary part are vanishing, the PT— symmetry of the system is spoilt and 
one obtains, in general, a spectrum not containing regions of purely real eigenvalues. 

For a Hermitian operator the spectrum is purely real. Possible level crossing 
(degeneration) points will be so-called diabolical points (DPs) [57], which are 
connected with diagonalizable spectral decomposition, where algebraic multiplicity 
Ua and geometric multiplicity Ug [11] of the degenerate eigenvalues A coincidqj, 
na{^) — nglX) and one finds a symmetry enhancement^. In addition to these diabolic 
degeneration points, for a non-Hermitian operator the occurrence of exceptional points 
(EPs) is possible and even generic. EPs are parameter configurations at which for a 
corresponding degenerate eigenvalue A the algebraic multiplicity na{X) is exceeding the 
geometric multiplicity, rig (A) < no (A). This is connected with a non-diagonalizable 
spectral decomposition of the operator (matrix) [10, 11], i.e., the formation of non- 
trivial Jordan-block structures [39, 59] and Jordan chains of algebraic eigenvectors 



^) Eigenvalues of diagonalizable matrices are called semi-simple (see e.g. [58]) and in case of 
^g{^) — "a(-^) = 1 simple. 

For a fc— fold DP a U (k) rotation symmetry occurs within the span of the degenerate eigenvectors. 
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(associated vectors) . Subsequently, we use the term mth-order EP for an EP which is 
associated with an mth-order Jordan block in the spectral decomposition. 

The EPs and DPs live on certain hypersurfaces Vj in the underlying parameter 
space M D Vj of the model. They are, in general, of various co-dimensions and form 
a so-called stratified manifold V = [jjVj (see, e.g., [60]). Depending on a concrete 
parameter perturbation the system may move along the stratified manifold V passing 
from one degeneration type to another one, or, more generically, it may escape from 
V so that the degeneration disappears and an EP or DP unfolds into non-degenerate 
eigenvalues. 

In the following sections we will analyze the occurrence and the unfolding of EPs 
for our Bose-Hubbard model. We will find that in the limit of vanishing particle 
interaction the model can be solved analytically and that there exist only two EPs of 
order A'^ -I- 1. In dependence on the direction of the perturbation in parameter space 
each of these EPs unfolds either into [{N + l)/2] eigenvalue pair^ (and, in case of 
N + 1 odd, into an additional zero-eigenvalue) or into eigenvalue triplets (third-order 
eigenvalue rings) and {N + 1) mod 3 single eigenvalues. 



3. The limit of vanishing interaction 

For vanishing interaction, c — 0, the Hamiltonian H in (jl2p is a complex linear 
combination of su(2) Lie algebra elements 

H = 2{-ijL^ + vL^) e sl{2,C). (19) 

Figure [1] shows real and imaginary parts of the eigenvalues of if as a function of 7 
for an example with = 11 particles and fixed v — 1. All eigenvalues are purely real 
for I7I < 1 and purely imaginary for I7I > 1. For I7I = 1 we observe a degeneracy 
of all eigenvalues, which will turn out to correspond to a full Jordan block, resp. 
an EP of order N-l-1. Furthermore the real and imaginary parts are axis symmetric 
with regard to both axes En = and r„ = — a behavior which is due to the high 
symmetry of the Hamiltonian (fT9|) . The eigenvalues can be easily obtained analytically 
by diagonalizing H in case of I7I ^ \v\ and bringing it to its non-trivial Jordan block 
form in case of I7I = \v\. 

For this purpose we make use of the Baker-Campbell-Hausdorff formula and the 
su{2) commutation relations (O to obtain the well known rotations and boosts over 
the algebra su{2) specifically needed for our analysis 

^-tOLyj^^^tOLy ^ coa{e)L,_+sm{e)L^, (20) 
e°'^yL^e~°'^^ = cosh(a)L3 +isinh(a)L^, (21) 
e-'^^-L^e'^^- = cos{e)Ly + sini9)L, . (22) 
With the help of the boost (PT|) and the identification 

cosh(a) = , 7 „ , sinh(a) = , ^ o (23) 
•\/7^ — ^7^ — 



the Hamiltonian p9)l can be reshaped as 

H = ~ 2i\/ 7^ — [cosh(a)L2 -I- i sinh(Q!)L2;] 

= -2^v'^ --i'^e'^^yL^e-~°'^y . (24) 



^) The notation [a] stands for the floor function which yields the highest integer less or equal a G R. 
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From the fact that the irrep ([9]) remains vahd for any complex extension of 5^(2) it 
follows the completeness of the basis vectors 

1= \l,m){l,m\ (25) 



7n— — l 



and with it 

{l,j\H\l,k) = -2Vv^-7^ {l,j\e"'''\l,m){l,m\LS,m'){l,m'\e~'''^y\l,k) 



I 



m.m' — — I 
I 



- 2Vv' - Y S-^^Srak (26) 



m— - 



where 

5„,fc:=(/,m|e-"^»|/,fc). (27) 

For I7I 7^ the boost functions remain finite so that \a\ < 00 and S 

remains regular. Therefore, S acts as similarity transformation which diagonalizes 
the Hamiltonian (fT9| . The eigenvalues of H can be read off from ([26]) as 



Xn^Er,- iVr, = U^V^ - (28) 
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where n = —N, —N + 2, . . . , TV — 2, iV. In the Uniit 7 — > ±v the diagonaUzation breaks 
down and the boost becomes singular: \a\ — > 00. Instead of eq. (PT|) we may use a 
rotation with 9 = -7r/2, i.e. 



(29) 



This allows us to represent the Hamiltonian as 

= 2ve'''^-/^L±e-"^''^^, (30) 

so that with Rnk {I, n|e^"-^^/^|/, k) we have 

(/,j|i/|?,fc) -2«i?-„i(/,m|L±|?,n)i?„fc, det(i?) ^ 0. (31) 

According to ([7]) only the elements {I, fc±l|L±|/, k) are non- vanishing. This means that 
the matrix {l,m\L±\l,n) is subdiagonal for and superdiagonal for L+. Matrices 
with all elements {(in}n=i their sub- or superdiagonals non-vanishing and all other 
elements equal to zero are similar to a single Jordan block Jjv -1-1 (0) with eigenvalue 
A = 0: 



/ 











0-2 







on-i 







a-N 




Q-ljAr+i(0)g, 



Q = diag ^l,ai,aia2, • • • , J^^flfc^ . (32) 

This is also immediately evident from the angular momentum operators: The ladder 
operators L± should only have a single eigenstate, namely "spin- up" resp. "spin- 
down", belonging to the eigenvalue A = 0. Since for I7I = |?;| the Hamiltonian is 
equivalent to L±, with total angular momentum N/2, these configurations correspond 
to a degenerate eigenvalue A = which is an EP of order iV-|- 1. Under variation of 7 it 
unfolds into [{N + l)/2] eigenvalue pairs (and in case of iV -|- 1 odd, into an additional 
zero-eigenvalue) according to P5|) . The result is illustrated in figured] We note that 
although such a pairwise unfolding of an (A^-l-l)th-order EP according to a square root 
law •^t)2 _ ^2 -y^rith different scaling pre-factors n appears physically rather generic, it 
is mathematically very special. Typically an (Af-l-l)th-order EP at A = unfolds under 
a smah perturbation |e| < 1 into an {N + l)-ring A w e'"^e^/(^+^\ k = 0, . . . , N 
of non-degenerate eigenvalues. This follows simply from an effective equation of the 
type A^"*"^ « e. Intuitively, one might expect the spectral behavior (|28|) with its 
decoupled square roots to result from a decomposition of H into [-^^^j^] second-order 
Jordan blocks rather than from a single {N + l)th-order Jordan block. The subject 
of the remaining sections will be to further clarify this specific (mathematically non- 
standard) spectral behavior of H. 

In contrast to the simple analytical structure of the model and its complete 
solvability in case of vanishing interaction c, the situation changes drastically 
for non-vanishing interaction. In this case the spectrum of the non-Hermitian 
PT— symmetric Hamiltonian (|12p is most efficiently studied with the help of numerical 
and perturbational techniques. In the next section we will present some numerical 
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results for non-vanishing interaction. Qualitative aspects of these findings are 
explained analytically with the help of perturbational techniques in section [51 

4. Numerical results for non-vanishing interaction 




-2-1012 -2-1012 



7 y 

Figure 2. Real- and imaginary parts of the eigenvalues An = E„ — iT„ of the 
Bose-Hubbard Hamiltonian l|12| l as a function of the non-Hermiticity 7 for d = 1, 
c = 0.1/N and N = 11 particles. 

Figures [2] and [3] demonstrate the typical spectral behavior of a Hamiltonian 
(|12p in dependence on the non-Hermiticity parameter 7 for fixed weak and moderate 
interaction strength c. In the concrete example, the twelve eigenvalues of an iV = 11 
particle system are shown. 

First we note that the non-vanishing interaction c reduces the symmetry of the 
spectrum: The symmetries with regard to sign changes of the complex coupling 
7 ^ —7 and of the imaginary spectral components r„ ^ — r„ are not altered. The 
first symmetry results from the fact that, regardless of the interaction, for symmetric 
modes it does not matter which of them is coupled to the source and which to the sink. 
Due to this 7 ^ —7 symmetry we can restrict our analysis to the parameter region 
7 > 0. The second symmetry (with regard to r„ ^ — r„) is a direct implication of the 
Krein space symmetry (|13p of H which causes non-real eigenvalues to occur always in 
complex conjugate pairs. The additional symmetry En ^ —En present for c = is 
lost in case c 0. This is due to the fact that for c 7^ the square root branch points 
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Figure 3. Real- and imaginary parts of the eigenvalues A,i = E„ — iTn of the 
Bose-Hubbard Hamiltonian II12I I as a function of the non-Hermiticity 7 for v = 1, 
c = 0.5/N and N = 11 particles. 




n 



Figure 4. Trajectories of the complex energy eigenvalues An = En — iTn of the 
Bose-Hubbard Hamiltonian I I12I I as a function of c with < cN < 0.1 (left) and 
< cN < 1 (right), for ■y = v, A'' = 11 particles. 
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(EPs) become shifted in different ways and do no longer coalesce in the parameter 
space = K'^ 9 (7,11,0). For |c| ^ 7^ ~ v'^ this behavior can be roughly 

described as an effective deformation of the spectral branches (1^51) of the type 

A„(c) = a„(c) + b„{c)n^v^ - 7^ - d„(c) a„(0) = d„(0) = 0, 6„(0) = 1 (33) 

which for c 7^ leads to 

E„ = Re [A„(c)] ^ -E^n = - Rc [A-„(c)] . (34) 

From figures [2] and [3] one clearly sees that for fixed v — 1 each of the two EPs of 
order twelve present in figure [T] at 7 = ±v splits up into six second-order EPs with 
different positions 7„ — ±\/w^ — dn(c). The fact that these special points are EPs is 
obvious from the graphics. One clearly sees, firstly, that these points are associated 
with transitions from real spectral branches to complex conjugate ones and that they 
are therefore branch points. Secondly, one observes that the lines which branch off 
from these points scale faster than linearly so that the points cannot be diabolical 
points which are connected with a linear scaling (see e.g. [22]). 

For increasing |c| the deviations from the spectrum at c = also increase (compare 
figures n [2] and©. 

Subsequently, it is convenient to relabel the absolute values of the EP positions 
|7n| =: 7fc=o-(n) in increasing order as 

< 71 < 72 < • • • • (35) 

Furthermore, we dub the {N + l)th-order EPs as mother EPs. 

Complementary information about the unfolding of these mother EPs for non- 
vanishing interaction strength c ^ can be gained by considering the behavior of the 
spectral branches at the former EP positions I7I — \v\. For this purpose the trajectories 
of the twelve eigenvalues in the complex plane have been plotted for interaction 
strengths c varying in the intervals < c < Cmax = 0.1/A^ and < c < c,„ax = 1/-^^ (see 
figure [4]). For definiteness we have chosen u = 1 so that the mother EP is localized at 
7 = u = 1. Obviously, for small |c| ^ \v\/N the twelve eigenvalues form three groups 
where four eigenvalues behave qualitatively almost identical — moving along one of 
the three lines in the directions ~ e~"^, e~''^^'"^ . This regular circle division with 
lines enclosing angles of 27r/3 in the complex plane clearly indicates on an unfolding 
of the type 

(-/j)'/V^ci/3, fc = 0,1,2, /, eM+,i = 1,2,3,4 (36) 

where three eigenvalues corresponding to the same fj can be understood as a triplet. 
For larger values of c (see , e.g., the right side of figure[3]) one eigenvalue of each triplet 
stays on the negative real axis, while the other two depart from their initial directions 
symmetrically with regard to the real axis as complex conjugate pairs. 

The triplet splitting is less obvious from the spectral branches depicted in figures 
[2]and[3l Nevertheless, these plots also clearly show that at the position -f = v — 1 four 
purely real eigenvalues and four complex conjugate pairs of eigenvalues are present. 
Obviously, the triplet splitting is closely connected with the fact that for small c ^ 
two of the second-order branch points (EPs) move to values 75 , 76 > 1 yielding in 
this way four real eigenvalues at 7 = 1, whereas the other four second-order EPs with 
7i, . . . ,74 move to positions < 7^ < 1, what results in the four complex conjugate 
pairs at 7 = 1. 

Further numerical investigation shows that the triplet unfolding of the mother 
EP for c ^ and I7I = \v\ is a generic feature of the system (|T2|) for arbitrary particle 
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number. In general, we find a splitting into ['^^j-'-] eigenvalue triplets and (iV+1) mod 3 
single eigenvalues. In section [51 this unfolding behavior will be analytically explained 
with the help of perturbational techniques. 




-2-1012 -2-1012 




-2-1012 -2-1012 



Y 7 

Figure 5. Real- and imaginary parts of the eigenvalues An = E„ — iTn of the 
Bose-Hubbard Hamiltonian II12I I as a function of the non-Hermiticity 7 for v = 1, 
A'^ = 11 particles and c = 1/A'^ (upper figures) resp. c = 2/N (lower figures). 




Figure 6. Real- and imaginary parts of the eigenvalues \„ = En — iT„ of the 
Bose-Hubbard Hamiltonian II12I I as a function of the non-Hermiticity 7 for v = 1, 
N = 11 particles and c = 80/Af. 



After getting a rough qualitative understanding of the unfolding of the mother 
EPs at I7I = \v\, c = we turn now to the system behavior for intermediate 



VT— symmetric Bose-Hubbard model 



14 




20 40 60 80 2 4 6 8 10 

cN cN 



Figure 7. (Color online) EP positions in the (c, 7)— half-plane 7 > for fixed 
V = 1 and N = 11. The figure on the right is a magnification of the small cN 
region. The EP positions are numerically approximated as the smallest positive 
values of 7 at which two formally degenerate resonance widths Fj j, differ by more 
than IQ-"'. 



and large interaction strengths c. Figures [5] and [B] give an impression about the 
corresponding spectral behavior and figure [7| about the associated location of the EPs 
in the (c, 7)— half-plane. Obviously, the movement of the second-order EPs remains 
qualitatively the same for intermediate values of c: two of the six EPs connected with 
the mother EP at 7 = u = 1 move away from 7 = 1 to positions 7fc > 1, whereas the 
other four 7^ remain located within the interval (0, 1) 9 7fc and move to smaller 7. In 
the limit of c 00 not only these four EPs tend asymptotically to the limiting value 
7 = (see figures [H] and [7|) , but also the EP 75 originally located at 7 > 1. The value 
7 = itself is only reached at c = 00, because for c < 00 and 7 = 6 the Hamiltonian 
H in ([8]) is a tridiagonal symmetric purely real matrix with non-vanishing elements 
on the superdiagonal. According to [61] (chapter 5, sections 36 and 37) such matrices 
have distinct (nondegenerate) eigenvalues. 

Figure O clearly shows the shrinking of real spectral "bubbles" which are defined 
by the intervals [—7^,7^], k < 5 for increasing c. The same fact is implicitly reflected 
in the EP behavior in the (c, 7)— half-plane (figure [7]). The 7— positions of the EPs 
and with them the sizes of the real spectral "bubbles" tend only asymptotically to 
zero. By zooming into figures [S] and [7] one would observe still remaining real-energy 
regions for all spectral branches shown in figure O We note that the smallest of those 
regions [—71,71] defines the sector of exact PT— symmetry (all eigenvalues are purely 
real). For c ^ oo this region together with the other intervals [— 7fc, 7fc], k < 5 shrinks 
to zero width (see figure [6]) . The shrinkage effect for purely real branches and the 
widening of branches with complex conjugate eigenvalues is similar to the well known 
effect of eigenvalue complexification for strong-coupling regimes in Sturm-Liouville 
type T^T— symmetric models [53,55]. 

In contrast to the EP accumulation 7^ ^ 0, fc < 5, for c ^ 00, the single EP pair 

±76 tends asymptotically to the limiting values ±7g°°'' = ±6. This is the manifestation 

of a generic result which holds for systems with odd particle numbers N and which is 

derived by perturbation techniques in section !??^ It states that for odd N and c 00 

the lowest two levels coalesce at EPs located at 7^^^! = ±v^^^^. In case of even N 

2 

all EP positions tend asymptotically to 7 = 0. 
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Let us restrict our attention (for symmetry reasons 7 ^ —7) to the region 
7 > 0. From the hmiting locations of the EPs for c — > and c — > 00 and the 
numerical results for c G [0, 00) we conclude that the second order EPs are located 
on two-dimensional hypersurfaces V„ = {(7,f,c) G A4 : j — 7„(?;,c)} in the three- 
dimensional parameter space Ai 3 {'-f,v,c). For c — these surfaces V„ coalesce at 
a common line Hn ~ lil^'^^'^) & -M : j = v,c = 0}. In the opposite limit of 
c ^ cxD all but one of the surfaces (in case of N odd) tend asymptotically to the plane 
Afa = {(7, ti, c) G : 7 = 0}. The remaining surface V w+i asymptotically approaches 

the plane A/'+ = {{j,v,c) & A4 : 7 = u^^^^}. Hence, the stratified manifold V C A4 
of EP degenerations for the present model comprises the two-dimensional surfaces Vn 
which intersect (coalesce) at the line p|„ V„ 



and its mirror images under the symmetry transformation 7 ^ —7. 

Summarizing the numerical findings we conclude that although the influence 
of the particle interaction might stabilize some of the energy levels, i.e. it might 
shift the occurrence of their imaginary parts to higher values of the non-Hermiticity 
parameter 7, the position of the first EPs monotonically decreases with increasing 
interaction strength c. Therefore the interaction always shrinks the region of unbroken 
PT symmetry. In the limit of infinitely strong interaction even an arbitrarily small 
complex perturbation destroys the reality of the spectrum. 

5. Perturbative results 

In this section we will investigate some of the features of the spectrum of 
the Hamiltonian (|12p for different limiting cases analytically using perturbational 
techniques. 

At first we will focus on the unfolding of the {N + l)th-order EP due to a 
weak particle interaction c. Since usual Rayleigh-Schrodinger perturbation theory 
breaks down for non-Hermitian Hamiltonians in the vicinity of EPs of the unperturbed 
operator, we resort to the Puiseux-Newton resp. Newton-polygon method [11,40,41], 
a perturbative technique for the roots of polynomials which works in the vicinity of 
degeneracies as well as around simple values. This method will allow us to analytically 
verify the triplet unfolding of the mother EP we found numerically in section [4l 

Finally, in the second part of this section, we will use standard Rayleigh- 
Schrodinger perturbation to understand the behavior of the spectrum in the strong- 
interaction regime. 

5.1. The limit of weak interaction 

Here we are going to analyze the unfolding of the {N + l)th-order EP when the 
interaction is switched on perturbatively, i.e., we consider the spectral behavior of the 
Hamiltonian IT^ at the EP 7 = u 




(37) 



n 



H = 2v{L^ - iL^) + 2cL\ 
for small interaction < |c| ^ bl/^- 



(38) 



VT— symmetric Bose-Hubbard model 



16 



As a first step, we S'C/(2)— rotate the Hamiltonian 
and = 7r/2 into the more convenient form 

H = 2v{L^ - iLy) + 2cLI 



with the help of eq. 



= 2wL_ 
=: 2wL_ 



L, — L(] 



Ln :— L+L- 



(39) 



In the particle number (angular momentum I) representation this Hamiltonian has a 
band diagonal structure of the type 



H 



( * 

* 

* 









where the second super- 
and the diagonal by Lq. 



\ 












(40) 



and subdiagonals are generated by the L\ and l?_ 



terms 

The perturbation matrix is an upper 3— Hessenberg matrix, 
i.e. a matrix with only zero entries below the three subdiagonals (including the main 
diagonal). Therefore the results of [62] apply, where the unfolding of the eigenvalues 
of Jordan blocks Jn(Ao) under perturbations by general upper fc— Hessenberg matrices 
has been analyzed. It has been shown that an nth-order EP typically splits into p 
rings of size fc and one of size r (if r 7^ 0), where n = pk + r, p = [^] and r = n modp. 
This means that for a small coupling parameter |c| ^ l^^l/-^ the EP will imfold as 



1 q=l, 



1. 



,P 
(41) 



The coefficients flg are specific model dependent constants whose moduli \aq\ define 
the scaling of the ring radii. The specific n = pk + r splitting behavior generalizes the 
well known results for the generic case {p = 1, r = 0) where the degenerate eigenvalue 
at the EP splits into a single ring of size n, i.e. where (in suitable reparametrization) 



A" 



holds, with the ring A = e 



J = 0, 



1, of size n as obvious 



solutions. Applying the results of [62] to the unfolding of our EP configuration we 
expect the formation of p — [^^3^] rings of size fc = 3 and possibly of a single ring 
of size r — 1 (single eigenvalue) or r = 2 (eigenvalue pair) depending on the concrete 
dimension N + 1 and r = (TV + 1) modp. 

Subsequently, we present an explicit derivation of this behavior which makes use 
of the specific structure of our model. The analysis will be based on the characteristic 
polynomial of the matrix Hamiltonian H (j39p as function of the interaction strength 
c, 



det(xi-H) = - ^ 



PM-k{c)X'' 



M N+1 



(42) 



fc=0 



which we will study with the help of the Newton-polygon technique [11,41]. This will 
allow us to derive the dominant fractional power fi of the unfolding A — uqc'^ + o{c'^) 
of the {N + l)th-order EP. 
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We start our analysis by noticing that the characteristic polynomial of any 
■PT— symmetric matrix Hamiltonian H has purely real coefficients. This follows 
straightforwardly from the fact that any T^T— symmetric operator can be represented 
as a purely real (possibly infinite dimensional) matrix [63]. 

The coefficients pM-fc(c) in flS]) can be recursively obtained with the help of the 
Le Verrier-Faddeev method [39] as 

1 * 

Pk = --^SjPk-j, fc = l,...,M, (43) 
i=i 

Sk ■= Tr(i/'=), po = -1 ■ 

For our purpose it is sufficient to extract the structure of the coefficients Pk{c) as 
polynomials in c. As basic input we first derive the corresponding traces Sfc(c) = 
Tr(iJ'^). These traces act in the angular momentum representation ([7]) or, equivalently, 
in the monomial representation ([9]). This means that only terms in contribute to 
Tr(_ff*^) which leave the angular momentum mode number m, and with it the monomial 
power, unchanged. Hence it will hold , e.g., Tr(L^) = Vfc G Z+, i.e. fc > 0, as well 
as Tr(L^Li) = Vfc 7^ j S Z+, but Tr(L^L^) ^ 0. In order to extract the power 
structure of Sfc(c) in c we may relate to L± their auxiliary commutative symbols 
L+ « ^, i_ K, ^^1. The terms can then be associated with the multinomials 

- E {\) v'-^e-'c^f: (-l)'e^(-') (44) 

(we omitted the irrelevant pre- factors 2 and —1/2 in front of v and c) where only the 
constant terms will contribute to the trace Tr(if'^). With the notation j :— i — l we 
find the terms as f° — ^^-^+"^3 , Combining the corresponding constraint k — l = 2j 
with the inequalities k>l>0,2l>i>0 one obtains j < I = k — 2j and, hence, the 
condition k > 3j or j < [|] . 

As result, we find the structure of the traces as 

m 

Skic,v) = Tt{H^) ^Y^afc'^-^'v'"'. (45) 

(k) 

The constant coefficients depend only on k and the particle number N (resp. the 
angular momentum I) and can be calculated with the help of general trace formulae 
for polynomials of angular momentum operators as discussed, e.g., in [64]. In our 
subsequent qualitative analysis we are only interested in the behavior of Sk{c,v) as 

(k) 

polynomial in c so that the concrete values of the non-vanishing coefficients aj are 
irrelevant. 

The trace formula (|45l) allows us to prove the following 
Theorem. The coefficients pk have the structure 
[-] 

Pk = i2df^^'~''^'' (46) 

(k) 

where d - are real constants which in general do not vanish. 
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Po = -1, Pi = -PqSi = si 
From and pS)) one finds pS)) by induction: 
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(47) 



fe=l 



n+1 [3j 



n 

E 

r=0 



fe=l i=0 



3 

E 

J=0 



^(fe) ^(n+l-fc) ^„+l-2(i+,) ^2(l+j) (48) 



^(n+l)^n+l-2r^2r- 



(49) 



5l" 



+1 



1 



n+1 



lEE E « 



(fc) ,(n+l-fc) 



<^r,i+j 



In passing from 



to 



fe=l (=0 j=0 

we used the fact that the total summation in goes 
over the three-dimensional discrete volume = {k,l,j : 1 < k < n + 1, < I < 
[f ] ) < J < [ "+1-*^ ]} and that Q is recovered by slicing it along fixed r — I + j and 
summing over the two-dimensional slices 



□ 



Summarizing the above results we conclude that for a given index fc it is sufficient 
to use the single index j as basic counting index, whereas the power of c in the 
polynomial terms is defined by the derived value I according to the relation: 

"^^ (50) 



l^k~2j 



J 



< 



As an illustration we list the c— dependence of the first six non-trivial coefficients 
Pk, fc = 1,...,6 of the characteristic polynomial Xij(A, c). According to (|46|) these 
coefficients contain terms of the following c'v^-'— monomial types 



fc = 1 


J=0, 


1, 




fc = 2 




2, 




fc = 3 


J=0, 1 = 


3, 








1, 




fc = 4 




4, 






3 = 1, 1 = 


2, 




fc = 5 


J=0, 1 = 


5, 






J = l, 1 = 


3, 




fc = 6 


J=0, 1 = 


6, 






J = 1, 1 = 


4, 






J=2, 1 = 


2, 


cV. (51) 


With the concrete monomial terms 


m 


at hand, it is now an easy task to obtain 


leading (dominating) exponent fix in the 


power series expansion 


A(c) = 






), 0<Aii </i2 < e, y^O (52) 
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for small |c| <C \v\/N. The method that we use is known as Puiseux-Newton diagram 
technique or Newton-polygon technique [11,40,41] and can be summarized as follows. 
One starts by substituting the ansatz (|52p with still unknown exponents /i^ and 
coefficients Ci into the characteristic polynomial 

X^[A(c),c] - A*^-pi(c)A*^-i pm-i{c)X-pm{c)^Q 

= (eic^- + • • -f - ifM-i<^"-' +■■■) (eic^^ + • • - f-' - 

ific"' +■■■) (eic^^ +•••)- ifoc"" +•••) = 0. (53) 

The dots in the various terms denote contribution of higher powers in c. For later 
convenience the lowest c— powers in the coefficients pk are numbered in reversed order 
as 

Pk{c) = fM-kc'""-" + o(c"^-^-'=) . (54) 

Since Xj^[A(c),c] as polynomial in c has to vanish, Xij[A(c),c] — 0, it should contain 
at least two terms in each power in c so that these terms can compensate one another. 
Taking into account that terms in the lowest power in c are the most dominating 
ones, one has to search for these lowest-power-terms and to fix the still undefined fj,i 
in such a way that these terms compensate. After fixing the minimal /ii one repeats 
the process for the next greater fj,2 and so on. In this way one can iteratively obtain 
the series expansion (15^ up to any required precision. The validity of the perturbative 
results is of course limited by the finite convergence radius of the perturbation series 
(EH). 

Returning to (|53p . one sees that fii as minimal exponent should be defined from 
the monomials 

efc^'^\ /A./-ief-^c'^"-^+(^'-^)^\ /M-2ef-2c'^^'-^+(''-'^^S...,/oc"'' (55) 

by pairwise identifying their powers 

flfc + fc^i.fc = aj + 3l^j,k, j 7^ k, j,k = 0,...,M. (56) 

In order to single out the relevant values of /ii one associates to each power ak + kfij,k 
a point Ak = (fc, a^) in the (fc, afc)— plane so that 

f^j,k = (57) 

K — J 

is just the sign- inverted slope of the line connecting the points Aj and Ak- As shown, 
e.g., in [11,40,41], the possible values of /xi can then be identified with those nj^k whose 
lines form the lower boundary of the convex hull of the points Ak- The corresponding 
graph is the so-called Newton polygon. Examples are depicted in figure El Points 
above this lower boundary of the convex hull will only contribute to higher order 
approximations . 

The coefficients ei in the series expansion (|52p can be obtained from reduced 
polynomials which are built from those leading-order terms of the characteristic 
polynomial ((53)) which correspond to points Ak located on the same lines of the 
Newton-polygon. In general, the lines comprise more than two points as it is visible, 
e.g., in figure m Here for TV = 5 the {^i = 1/3)— line comprises three points and for 
= 10 the (/ii = 1)— line also comprises three points and the (/ii = 1/3)— line four 
points. 

Applying the described Puiseux-Newton technique to our concrete characteristic 
polynomial we read off from eq. (|46p that the lowest c— powers in the coefficients Pk{c) 
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Figure 8. Puiseux-Newton-diagram for N = 5 particles (left) and N = 10 
particles (right). 



have the form d\'l\c''~^m With Af = TV + 1 and the reverse numbering (|54p we find 
the relevant points Ak in the (fc, afe)— plane as 

"fcl 



AN+i-k ^ [N+l-k,k-2 



k = 0,...,N + 1. (58) 



In case of iV = 5 this gives, for example, the seven points 

^ = (0,2), v4i = (l,3), A2 = (2,2), ^3^(3,1) 

^4 = (4,2), A = (5,1), A6 = (6,0) (59) 

which are depicted in the left graphics of figure [H Both of the two graphics in this 
figure show a typical modulo-threc-ratchet-structure of the points Ak- The lower 
boundary of the convex hull of these points is always formed by one (long) line of 
slope —1/3 (and corresponding fii = 1/3) which connects ["^^^-^-j + 1 points Ak and 
possibly a second (short) line of slope —1 (and = 1) which connects the first two 
or three leftmost points. We arrive at the result that in our model only dominant 
eigenvalue scalings of the type 

Acxcl/^ Aocc (60) 

are possible. Obviously, these eigenvalues will form rings of size three (triplets) and of 
size one (single eigenvalues). The specific size-one rings can be considered as atypical 
cases in the scheme of [62] . They can be attributed to the specific substructure of the 
Hessenberg type perturbation matrix. Moreover these X c terms can be identified 
as higher-order corrections — what is clearly visible by setting c^/'^ =: e so that 
c = e^. This is a specific output of the Newton-diagram technique which yields 
the dominant terms for all roots of a polynomial. In our concrete case this means 
that the 3— rings correspond to dominant lowest-order scaling c^/^, whereas the single 
eigenvalues remain unperturbed in this lowest order. Their dominant terms start only 
with the higher-order A ^ c contributions. 

Let us now illustrate the general theoretical results by more explicit calculations 
for models with N = 5 and TV = 10 particles. For the N = 5 model the corresponding 
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Hamiltonian has the form 



H 
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5i7 + 
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(61) 



In case of 7 = f the characteristic polynomial reads 
= det(iJ- A/) 

1743 



= A*^ - 35cA^ 



/82831 
V 16 
50625 



4645 



58275 
IT' 



+ 27280z;2c3 



A 



, . 30600t;2c4 + 6400i;V. (62) 

Obviously, A = is the only root for c — 0. For |c| > the coordinates of the points 
(|59p in the Puiseux-Newton diagram are given by the exponent of A and the minimal 
exponent of c of each summand. Under the conditions described above one finds one 
straight line with the slope —1/3 and associated /ii = 1/3. The three summands 
corresponding to the three points on this straight line must compensate each other 



6400c2w4 + 448^;^cA3 + A^ = 0, A = eic^/^ 
This yields two solution triplets for the coefficient ei: 
e\^'> = -y/U.77vi, -^433. 23?; ^ 



(63) 



(2) 

ef'> = \/14.77v 



3e 



(5) 

J6) 



v/433.23?; 



3e 3 . 



(64) 



Figure [S] shows the real parts (left graphics) and the imaginary parts (right 
graphics) of the numerically evaluated eigenvalues of the Hamiltonian for a small 
interaction strength c = 0.1/7V. On the line 7 = we find two pairs of complex 
conjugate eigenvalues with positive real parts and two purely real negative eigenvalues, 
in perfect agreement with the qualitatively predictions of the first order perturbation 
coefficients ([M)) . 

The numerically obtained eigenvalue trajectories in the complex plane for c G 
(0, 1) and fixed j = v are shown in figure (TUl For small c ^ we find again the 
typical 3— rings. Their radii are given by the two different absolute values of the 
coefficients (|64|) . For larger values of c the cubic-root-behavior becomes deformed 
by higher order corrections what is clearly visible as a deviation from the straight 
lines ~ e~"^^'"r. We note that the exact trajectories can be approximated in terms 
of a series expansion with arbitrary precision. Such series expansions are well defined 
over parameter regions which lie within the convergence radius of the series and break 
down when the next located algebraic singularity (branch point, EP) is reached. (The 
distance to the next EP defines the convergence radius (see, e.g., [65]).) 

Let us now turn to the model with iV = 10 particles. The corresponding Puiseux- 
Newton-diagram in figure [5] shows again the typical modulo-three-ratchet-structure. 
The two straight lines which form the lower boundary of the convex hull of the point 
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Figure 9. Real- and imaginary parts of the eigenvalues X„ = E„ — iTn of the 
Bose-Hubbard Hamiltonian l|12| l as a function of the non-Hermiticity 7 for v = 1, 
N = 5 particles and c = 0.1/A'^. 




Figure 10. (Color online) Trajectories of the complex energy eigenvalues 
An = En — ir„ of the Bose-Hubbard Hamiltonian II12I I as a function of c with 
< cAf < 1 (red), < cA^ < 0.5 (yellow) and < cA^ < 0.1 (orange), for = v, 
N = 5 particles. 



set have slopes —1/3 and —1 so that the llth-order EP unfolds with dominant scaling 
powers ^ = 1/3 and = 1, i.e. as A ~ c^/"^ and A ^ c. Reinserting these solutions 
into the corresponding characteristic polynomial yields one ninth-order equation for 
the coefficient ei in A « Cic^^^. This ninth-order equation reduces to a cubic equation 
in ef and gives the three different values |ei| as scaling parameters (ring radii) of 
the three triplets which comprise the first nine coefficients e^^'', . . . ,e^^''. The result 
is similar to eq. (|64|) only with three triplets instead of two. The quadratic equation 
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y 

Figure 11. Real- and imaginary parts of the eigenvalues X„ = En — iVn of the 
Bose-Hubbard Hamiltonian H12| l as a function of the non-Hermiticity 7 for v = 1, 
N = 10 particles and c = 0.1/A'^. The lower figure shows a magnification of the 
real parts near the central EP. 




Figure 12. Trajectories of the complex energy eigenvalues A„ = En — iTn of the 
Bose-Hubbard Hamiltonian H12| l as a function of c with < cN < 0.1, for -y = v, 
N = 10 particles and a magnification of the innermost region (right). 
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for the remaining two coefRcients ,e\ in the linear scaling law A ~ eic is easily 
derived by computer algebra and takes the explicit form 

- 46423756800e? + 2410418995200ei - 33581039616000 = 0. (65) 

It yields the coefficients 

145304 //145304\^ 4048640 „^ „. 

Real and imaginary part of the spectrum as well as the unfolding of the llth-order EP 
at c = and 'y — v — I are shown in figures [TT] and [121 Clearly visible in figure [12] are 
the three 3— rings and the two linearly scaling single eigenvalues. Obviously, eigenvalue 
shifts induced by the linear scaling are much smaller than the shifts of the 3— rings. In 
figures [m and [T2l these higher-order (A ~ c) corrections are only visible in the zoomed 
graphics. In leading c^^^— order approximation, the EP related to the A ~ c branches 
remained fixed at the original position j = v — 1. 

Summarizing we conclude that for arbitrary particle number N the Newton- 
diagram at the (A^-l-l)th-order EP shows a modulo-three-ratchet-structure with regard 
to the unfolding due to increasing interaction strength c — like in figure [H] For non- 
vanishing c the EP unfolds into ['^^j-'-] eigenvalue triplets forming regular 3— rings in 
the complex plane with dominant scaling behavior of the type c^^^ for |c| ^ l"*^!/-^- 
The remaining {N + 1) mod 3 single eigenvalues depend linearly on c. Furthermore we 
find for small interaction strength c and an original {N + l)th-order EP at I7I = |u| 
that roughly 2/3 of the occurring second-order EPs are located in parameter regions 
< I7I < \v\ and roughly 1/3 in the region \v\ < \j\ what confirms the numerical 
results of section 

It remains to emphasize that the unfolding of (A^+l)th-order EPs into triplets (for 
7^ = w^) has its origin in the effective 3— Hessenberg form of the perturbation matrix 
(|39|) . (|40|) . This suggests to reinterpret the square-root spectrum of the exactly solvable 
(c = 0)— model from section[3]as EP unfolding under perturbation by a 2— Hessenberg 
matrix. Indeed, representing the Hamiltonian (fT9|) as 

H ^2v{L^,-iL^)~2il\L^ A 7 - v (67) 

and performing an SU (2) rotation as for (|39p it takes the structure 

H = 2vL_- A{L+- L_). (68) 

Due to the fact that i_ is of Jordan block type the Hessenberg perturbation theory 
of [62] is applicable. The perturbation matrix A(L+ — L_) has non- vanishing 
entries only on the first sub- and superdiagonals and it is therefore of 2— Hessenberg 
type. According to [62] the (iV -I- l)th-order EP at A = unfolds then under this 
2— Hessenberg type perturbation into [■^^^^] eigenvalue pairs and, for TV even, into one 
additional single eigenvalue. Obviously, this prediction is in complete agreement with 
the exact result ([25]) for the spectrum which shows a square root (pairwise) unfolding 
of the mother EP and an additional single eigenvalue A = in case of even N. 

The EP unfolding according to the Hessenberg perturbation type is straightfor- 
wardly extendable to Hamiltonians with "f — v and higher-order perturbations in Lz 

H = 2v{L^-iLz) + cLl k>2. (69) 

In this case an SU{2) rotation leads to 

H = 2vL- c{L+-L-f. (70) 
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Figure 13. Eigenvalue trajectories of a (A^ = 4)— particle Bose-Hubbard 
Hamiltonian 1(6911 with fixed w = 1 for fc = 3 and < c < 0.1 /N^ (left, 4-ring) 
and for fc = 4 and < c < 0.1/N^ (right, 5— ring). 



For 1 < fc < the perturbation matrix {L+ — L-Y is of (fc-l- 1)— Hessenberg type 
so that the (A^ + l)th-order EP wiU unfold into [-f^] eigenvalue rings of size fc-l- 1 and 
r = {N + 1) mod (fc -I- 1) eigenvalues which will be grouped in one or several smaller 
rings. Figure [13] shows the EP unfolding for {N — 4)— particle Hamiltonians (l69|) with 
fc = 3 and k — A, small c and fixed 7 = w = 1. Clearly visible are the 4— ring and 
5— ring eigenvalue structures. For fc > TV also the lowest left matrix element becomes, 
in general, nonvanishing Hn+i,i ^ so that for these fc the mother EP will unfold 
into a single {N + 1)— ring of eigenvalues. 

5. 2. The limit of strong interaction 

To understand the limit of strong interaction analytically, we can apply ordinary 
perturbation theory with 7 and v being the small parameters: 

H = 2cLl -2i-fL,^ + 2vL^^ . (71) 

Ho Hi 

In general the eigenstates of Hq are doubly degenerate, which can be seen in the 
standard basis ([7|) 

2cLl\l,±m^) =2cml\l,±m,), (72) 

with I — An exclusion is the eigenstate with = in a system with even particle 
number N. This state is not degenerate. 

In lowest-order approximation, the perturbed degenerate energy levels are given 

as 

Ec^Ea+Ei, (73) 
where the corrections Ei can be calculated from the perturbation matrix 

W = {l,m::\Hi\l,m'^) m,,m'^ = ±\m,\. (74) 
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We start with even particle numbers N. In this case the matrix W is diagonal and 
the energy corrections read 

El = -2ijm,. (75) 

Obviously, the state to^ — remains unperturbed in lowest order approximation. 

For odd particle numbers N the correction ([75]) holds as well, except in the case 
of \mz\ = 1/2 where the perturbation matrix is nondiagonal 

The eigenvalues of this matrix are 

Ei^±sjv^ (^^^-^2 (77) 

and we see that there occur two second-order EPs at 7 = i-^^iif. The corrections 
El for the states \mz\ = \ are purely real for I7I < and purely imaginary for 

\l\>\v\^. 

The exact spectrum for c = and the numerical studies for c 7^ show that 
pairwise complex conjugate eigenvalues of H occur for large values of I7I . In connection 
with the purely imaginary perturbative corrections Ei for states with \mz\ > 1/2 this 
implies that for large |c| 3> I'l'l/^ all EPs involving these states must have tended to 
7^0. 

Summarizing we conclude that in the limit c — > cx) there exists one zero-eigenvalue 
state with = for A'^ even and a pair of real eigenvalues for states jm^l = 1/2 in 
the parameter region I7I < \v\^ of a model with N odd. All remaining eigenvalues 
come as complex conjugate pairs. 

Therefore these perturbative results prove the numerical observations of section 
HI that in the limit of ultra-strong interaction all EPs of an A^— even model are located 
at 7 = 0, whereas in an TV— odd model two of the EPs can be found at 7 = ±.v^^^. 

6. Conclusion and Outlook 

We studied the spectrum of a non-Hermitian PT— symmetric two-mode Bose-Hubbard 
Hamiltonian, a system modeling an A^— particle Bose-Einstein condensate in a double 
well potential containing a sink in one of the wells and a source of equal strength in the 
other. While for vanishing particle interaction there exists only one pair of EPs of order 
A'^-l-l, the interplay of non-Hermiticity and particle interaction leads to a characteristic 
unfolding of these EPs into 3— rings of eigenvalues and the occurrence of a series of EPs 
of order two. This numerically observed scenario has been analytically understood 
using the Puiseux-Newton perturbation technique. Furthermore the case of strong 
particle interaction was described by ordinary Rayleigh-Schrodinger perturbation 
theory. 

Further investigations concerning, e.g., the positions of the EPs as well as their 
influence on the system dynamics remain tasks for future research. 

Another challenge is the investigation of the large N limit of the present model, 
resp. the so called mean-field approximation. In the Hermitian case this mean- 
field approximation is usually achieved by replacing the bosonic field operators by 
c-numbers, the condensate wave functions, yielding the nonlinear Schrodinger equation 
resp. Gross-Pitaevskii equation. This approach is closely related to a classicalization. 
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In a number of recent papers consequences of the classical nature of the mean-field 
approximation are discussed and semiclassical aspects are introduced [33,35,36,66-68]. 
For a two-mode system even the eigenenergies and eigenstates of the many particle 
system could be reconstructed approximately from the mean-field system in a 
semiclassical approximation with astonishing accuracy [69]. While there are some 
investigations concerning an heuristically introduced non-Hermitian generalization of 
discrete nonlinear Schrodinger equations [37,70-74] a careful derivation of a mean-field 
approximation starting from a non-Hermitian many particle system was lacking in the 
past and will be the subject of a separate paper [75]. 
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